# ____________________________________________________________
# Replication for Appendix Results
# ____________________________________________________________

rm(list = ls())
source("code/_preamble.R"); preamble()
source("code/_functions.R")

# load data

# election data
kompally_dta <- read_csv("data/election/kompally.csv")
ward_dta <- read_csv("data/election/kompally_ward.csv")
ge_dta <- read_csv("data/election/kompally_general_election.csv")
candidate_dta <- read_csv("data/election/kompally_candidate.csv")
# frt-use data
frt_worldwide <- read_csv("data/frt-use/frt_worldwide.csv")
frt_india <- read_csv("data/frt-use/frt_india.csv")
# census data
census <- read_csv("data/census/census_select_districts.csv")
# google-trends data
google_trends <- read_csv("data/google-trends/frt_google_trends.csv")
# maps
india_state_shp <- st_read("data/maps/india_state.shp")
te_shp <- st_read("data/maps/telangana.shp")

sims <- 100000
sims <- 10

# ____________________________________________________________
# Appendix 1: Context ----
# ____________________________________________________________

# Figure A1 --------------------------------------------------------------------
# Use of FRT Worldwide (A1.1) ---
world <- ne_countries(scale = "medium", returnclass = "sf")
world <- world %>% select(name, iso_a3, geometry) %>% filter(!is.na(iso_a3))

map <- left_join(world, frt_worldwide, by = c("iso_a3" = "code"))

map <- map %>% mutate(
  frt_status = factor(frt_status, levels = c("In use", "Approved for use (not implemented)", "Considering technology",
                                             "Banned", "No evidence of use"))
)

# calculate share of countries where FRT is in use
map %>% filter(frt_status == "In use") %>% nrow() / 
  map %>% filter(!is.na(frt_status)) %>% nrow()

cbp1 <- c("dodgerblue3", "deepskyblue1", "lightblue1", 
          "goldenrod3", "gold", "#D55E00", "#E69F00")

# plot
ggplot(data = map) +
  geom_sf(aes(fill = frt_status)) + 
  scale_fill_manual(values = cbp1, na.value = "white") +
  theme_map() + theme(legend.position="bottom",
                      legend.title=element_text(size=6),
                      legend.text=element_text(size=6)) + labs(fill = "FRT Status") 
ggsave("figures/world.pdf")


# Figure A2 --------------------------------------------------------------------
# Use of FRT in India (A1.2) ---

frt_india_map <- india_state_shp %>% 
  left_join(frt_india, c("ST_NM" = "state")) 

# calculations
frt_india %>% pull(frt_systems) %>% sum()
frt_india_map %>% filter(frt_systems > 0) %>% nrow()

cbp1 <- c("#999999", "#009E73", "#CC79A7", "#56B4E9", 
          "#F0E442", "#0072B2", "#D55E00", "#E69F00")

frt_india_map <- sf::st_simplify(frt_india_map, preserveTopology = TRUE, dTolerance = .01)

ggplot() +
  #geom_sf(data = india_country_map, color = "gray6", fill = NA) +
  #geom_sf(india_country_map, aes(fill = "white")) + 
  geom_sf(data = frt_india_map, aes(fill = frt_systems)) + 
  scale_fill_gradient(low = "white", high = "steelblue4") + 
  theme_map() + theme(legend.position="bottom",
                      legend.title=element_text(size=6),
                      legend.text=element_text(size=6)) + labs(fill = "Number of FRT Systems") 
ggsave("figures/indiafrt.pdf")

# Figure A3 --------------------------------------------------------------------
# Kompally Map (A1.3) ---

label <- cbind(c("Hyderabad", "Kompally"), c(78.48532998335422, 78.4848918485281), c(17.38823774025801, 17.537576520002215)) %>% 
  as_tibble() %>% 
  mutate(name = V1,
         long = as.numeric(V2),
         lat = as.numeric(V3)) %>% 
  select(name:lat)

te_shp <- sf::st_simplify(te_shp, preserveTopology = TRUE, dTolerance = .01)

ggplot() +
  geom_sf(data = te_shp, fill = "azure2", color = "gray62") + 
  geom_point(data = label, aes(x = long, y = lat)) +
  geom_text_repel(data = label, aes(x = long, y = lat, label = name), size = 5,
                  fontface = "bold", segment.size = 0.5, nudge_x = c(-0.85, 0.75), nudge_y = c(-0.5, 0.5)) + 
  theme_map()
ggsave("figures/kompally.pdf")

# Figure A4 --------------------------------------------------------------------
# Compare District with Kompally to Others (A1.3) ---

distrist_list <- c("083", "086", "076", "088", "139",
                   "140", "141", "518", "537", "584",
                   "583", "631", "481", "475", "471",
                   "473", "483", "482")

# Get urban areas
urban <- census %>% 
  filter(Level == "DISTRICT") %>% 
  filter(TRU == "Urban") %>% 
  select(District, urban_pop = TOT_P)

# Muslim population per district obtained from: https://www.censusindia.co.in/districts 
muslim_pop <- data.frame(District = distrist_list,
                         muslim_share = c(0.0086, 0.0468, 0.0311, 0.0893, 0.2798,
                            0.2535, 0.1308, 0.1919, 0.1166, 0.1056,
                            0.0931, 0.0613, 0.0739, 0.0625, 0.0670,
                            0.0412, 0.1150, 0.1199)) %>% 
  
  as_tibble()

# create district level variables
districts <- census %>% 
  filter(District %in% distrist_list) %>% 
  filter(Level == "DISTRICT") %>% 
  filter(TRU == "Total") %>% 
  left_join(urban) %>% 
  left_join(muslim_pop) %>% 
  mutate(sex_ratio = TOT_M / TOT_F,
         sc_share = P_SC / TOT_P,
         st_share = P_ST / TOT_P,
         urban = urban_pop / TOT_P,
         male_literacy = M_LIT / TOT_M,
         female_literacy = F_LIT / TOT_F,
         workers_share = TOT_WORK_P / TOT_P,
         agr_share = MAIN_AL_P / TOT_P,
         group = ifelse(Name == "Rangareddy", "Ranga Reddy", "Other Districts Near Major Cities")) %>% 
  select(Name, group, muslim_share, sex_ratio:agr_share)

# list of variables
var_list <- districts %>% 
  select(muslim_share:agr_share) %>% 
  names()

# for loop to generate summary statistics
mylist <- list()
for (i in 1:length(var_list)){
  
  t <- summarySE(districts %>% filter(!is.na(group)),
                 measurevar = var_list[i], groupvars = c("group"), na.rm = T) %>% 
    mutate(var_list = var_list[i])
  
  mylist[[i]] <- t
}
# Note: Ranga Reddy will not produce any CIs since it's just one observation

# Prepare plot data and variable names
plot_dta <- do.call("rbind", mylist) %>% as_tibble() %>% 
  mutate(group = factor(group, levels = c("Ranga Reddy", "Other Districts Near Major Cities")))

vl <- data.frame(var_list = var_list, 
                 var_label = c("Percent Muslim", "Sex Ratio", "Percent SC", "Percent ST", "Percent Urban",
                               "Male Literacy", "Female Literacy", "Percent of Workers", "Percent in Agriculture"
                 )) %>% as_tibble()

plot_dta <- plot_dta %>%
  mutate(value = ifelse(var_list == "sex_ratio", value, round(value*100, 1)),
         ci = ifelse(var_list == "sex_ratio", ci, ci * 100)) %>% 
  left_join(vl)

# Plotting
ggplot(plot_dta, 
       aes(x = group, y = value)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin = value - ci, ymax = value + ci),
                width=.2,                    
                position=position_dodge(.9)) +
  facet_wrap(~var_label, scales = "free_y") + 
  geom_text(aes(label = round(value, 1)), size = 4, vjust = 2, family = "mono", fontface = "bold", color = "white") + 
  xlab("") + ylab("") + 
  theme_bw() 
ggsave("figures/district_comparison.png", height = 10, width = 14)

# ____________________________________________________________
# Appendix 2: Data and Identification ----
# ____________________________________________________________

# Table A2 --------------------------------------------------------------------
# Summary Statistics (A2.1) ---

# select variables
ss <- kompally_dta %>% 
  select(turnout, facial_recognition, male_turnout, female_turnout,
         total_eligible, male_eligible, female_eligible, pfemale, n_party,
         percentmuslim_method1, percentmuslim_method2, women_res, stsc_res, bc_res, turnout_ge)

# create summary stats table
stat.desc(ss) %>% 
  round(digits = 2) %>% 
  t() %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column("variable") %>% 
  mutate(var_label = c("Turnout", "Facial Recognition", "Male Turnout", "Female Turnout",
                       "Total Eligible", "Male Eligible", "Female Eligible", "Share of Female", "Number of Parties",
                       "Percent Muslim (M1)", "Percent Muslim (M2)", "Gender Reservation", "SC / ST Reservation",
                       "BC Reservation", "General Election Turnout")) %>% 
  select("Variable" = var_label, "Mean" = mean, "SD" = std.dev, "Min" = min, "Max" = max, "N" = nbr.val) %>% 
  xtable() # %>% 
  # Kable() %>% 
  # KableSave(fname = "summarystat")

# Figure A5  --------------------------------------------------------------------
# Plot DV (Turnout) (A2.1) ---
ggplot(kompally_dta, aes(x = turnout)) + 
  geom_histogram(aes(y=..count..), bins = 12, colour="black", fill="gray58")+
  geom_vline(aes(xintercept=mean(turnout)),
             color="blue", linetype="dashed", size=1) +
  theme_bw() + xlab("Turnout") + ylab("Count")
ggsave("figures/turnout.png")

# Figure A6  --------------------------------------------------------------------
# Citizen's Knowledge of FRT in Google Trends (A2.2) ---

google_trends <- google_trends %>% 
  mutate(date = mdy(date),
         week = floor_date(date, unit = "week")) %>% 
  group_by(week) %>% 
  dplyr::summarise(frt = mean(frt)) %>% 
  filter(week > as.Date("2020-01-01")) %>% 
  filter(week < as.Date("2020-02-07")) 

ggplot(google_trends) + 
  geom_line(aes(x = week, y = frt), lwd = .75) +
  geom_point(aes(x = week, y = frt), lwd = 2, color = "blue") +
  scale_x_date(breaks=c(as.Date("2020-01-05"), as.Date("2020-01-12") , as.Date("2020-01-19"), 
                        as.Date("2020-01-26"), as.Date("2020-02-02"))) + 
  geom_vline(xintercept = as.Date("2020-01-18"), linetype="dotted") + 
  geom_vline(xintercept = as.Date("2020-01-22"), linetype="dotted") +
  geom_vline(xintercept = as.Date("2020-01-25"), linetype="dotted") +
  theme_classic(14) +
  ylab("Average Google Trends Index of Interest") + 
  xlab("Date") + 
  ylim(c(0,20)) +
  geom_text(aes(x = as.Date("2020-01-18"), y = 17, label = "TSEC Circular"), color="blue", size = 4) + 
  geom_text(aes(x = as.Date("2020-01-22"), y = 15, label = "Election Day"), color="blue", size = 4) + 
  geom_text(aes(x = as.Date("2020-01-25"), y = 18, label = "Results"), color="blue", size = 4)
ggsave("figures/frt_googletrends.png")

# ____________________________________________________________
# Appendix 3: Additional Analysis and Robustness ----
# ____________________________________________________________

# Table A3  --------------------------------------------------------------------
# Effect of Gender Turnout (A3.2) ---

m1 <- felm(male_turnout ~ facial_recognition | 0 | 0 | 0, 
           data = kompally_dta)
m1 %>% summary(., robust = T)

m2 <- felm(female_turnout ~ facial_recognition | 0 | 0 | 0, 
           data = kompally_dta)
m2 %>% summary(., robust = T)

# Randomization Inference to get exact p-values
set.seed(1234)
m1_p <- conduct_ri(
  formula = male_turnout ~ facial_recognition,
  declaration = declare_ra(N = 36, m = 10),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = kompally_dta,
  sims = sims
)
m1_p

m2_p <- conduct_ri(
  formula = female_turnout ~ facial_recognition,
  declaration = declare_ra(N = 36, m = 10),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = kompally_dta,
  sims = sims
)
m2_p

m1_pval <- summary(m1_p)$two_tailed_p_value %>% round(2)
m2_pval <- summary(m2_p)$two_tailed_p_value %>% round(2)

stargazer(m1, m2,
          keep = c("facial_recognition"),
          covariate.labels = c("Facial Recognition Technology"),
          se = list(m1$rse ,m2$rse),
          star.cutoffs = NA,
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          report=('vc*s*p'),
          dep.var.labels   = c("Male Turnout", "Female Turnout"),
          keep.stat = c("n"),
          title = "Effect of Facial Recognition Software on Voter Turnout",
          add.lines = list(c("Mean of DV", round(mean(kompally_dta$male_turnout), 2), round(mean(kompally_dta$female_turnout), 2))),
          float = F,
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:gender-effect",
          out = "tables/gender-effect.tex") 

# Replace regular p-values with exact p-values (Fisher)
paste0("& [", m1_pval, "]", "$^{}$", " & [", m2_pval, "]", "$^{*}$", "\\")


# Table A4  --------------------------------------------------------------------
# Controlling for General Election Turnout (A3.2) ---

m1 <- felm(turnout ~ facial_recognition + turnout_ge | 0 | 0 | 0, 
           data = kompally_dta)
m1 %>% summary(., robust = T)

m2 <- felm(turnout ~ facial_recognition + pfemale + total_eligible + 
             n_party + any_res + percentmuslim_method2 + turnout_ge | 0 | 0 | 0, 
           data = kompally_dta)
m2 %>% summary(., robust = T)

d <- kompally_dta %>% filter(!is.na(turnout_ge)) 

set.seed(1234)
m1_p <- conduct_ri(
  formula = turnout ~ facial_recognition + turnout_ge,
  declaration = declare_ra(N = 26, m = 7),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = d,
  sims = sims
)
m1_p

m2_p <- conduct_ri(
  formula = turnout ~ facial_recognition + pfemale + total_eligible + 
    n_party + any_res + percentmuslim_method2 + turnout_ge,
  declaration = declare_ra(N = 26, m = 7),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = d,
  sims = sims
)
m2_p

m1_pval <- summary(m1_p)$two_tailed_p_value %>% round(2)
m2_pval <- summary(m2_p)$two_tailed_p_value %>% round(2)

stargazer(m1, m2, 
          keep = c("facial_recognition"),
          covariate.labels = c("Facial Recognition Technology"),
          se = list(m1$rse, m2$rse),
          star.cutoffs = NA,
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          report=('vc*s*p'),
          dep.var.labels.include = FALSE,
          keep.stat = c("n"),
          title = "Effect of Facial Recognition Software on Voter Turnout",
          add.lines = list(c("Mean of DV", rep(round(mean(d$turnout), 2), 2)),
                           c("Previous Turnout Control", "Y", "Y"),
                           c("Election Controls", "N", "Y"), 
                           c("Polling Station Controls", "N", "Y")),
          float = F,
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:prevturnout",
          out = "tables/prevturnout-effect.tex") 

# Replace regular p-values with exact p-values (Fisher)
paste0("& [", m1_pval, "]", "$^{*}$", " & [", m2_pval, "]", "$^{}$", "\\")


# Table A5  --------------------------------------------------------------------
# Controlling for Post-Treatment Variables (3.4) ---

# Get ENOP
enop <- candidate_dta %>% 
  group_by(ward) %>% 
  dplyr::summarise(enop = enp(percent_vote),
                   margin = sum(diff_win_second, na.rm = T))

# Models
kompally_dta <- left_join(kompally_dta, enop, by = c("ward" = "ward"))
rm(enop)

m1 <- felm(turnout ~ facial_recognition + enop + margin  | 0 | 0 | 0, 
           data = kompally_dta)
m1 %>% summary(., robust = T)

m2 <- felm(turnout ~ facial_recognition + pfemale + total_eligible + 
             n_party + any_res + percentmuslim_method2 + enop + margin | 0 | 0 | 0, 
           data = kompally_dta)
m2 %>% summary(., robust = T)

set.seed(1234)
m1_p <- conduct_ri(
  formula = turnout ~ facial_recognition + enop + margin,
  declaration = declare_ra(N = 36, m = 10),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = kompally_dta,
  sims = sims
)
m1_p

m2_p <- conduct_ri(
  formula = turnout ~ facial_recognition + pfemale + total_eligible + 
    n_party + any_res + percentmuslim_method2 + enop + margin,
  declaration = declare_ra(N = 36, m = 10),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = kompally_dta,
  sims = sims
)
m2_p

m1_pval <- summary(m1_p)$two_tailed_p_value %>% round(2)
m2_pval <- summary(m2_p)$two_tailed_p_value %>% round(2)

stargazer(m1, m2, 
          keep = c("facial_recognition"),
          covariate.labels = c("Facial Recognition Technology"),
          se = list(m1$rse, m2$rse),
          star.cutoffs = NA,
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          report=('vc*s*p'),
          dep.var.labels.include = FALSE,
          keep.stat = c("n"),
          title = "Effect of Facial Recognition Software on Voter Turnout",
          add.lines = list(c("Mean of DV", rep(round(mean(kompally_dta$turnout), 2), 2)),
                           c("Post Treatment Controls", "Y", "Y"),
                           c("Election Controls", "N", "Y"), 
                           c("Polling Station Controls", "N", "Y")),
          float = F,
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:posttreatment",
          out = "tables/posttreatment-effect.tex") 

# Replace regular p-values with exact p-values (Fisher)
paste0("& [", m1_pval, "]", "$^{*}$", " & [", m2_pval, "]", "$^{*}$", "\\")

# Table A6  --------------------------------------------------------------------
# Covariate Adjustment, Lin (2013) (3.5) ---
m1 <- lm_lin(turnout ~ facial_recognition,
             covariates = ~ pfemale + total_eligible + n_party + any_res + percentmuslim_method2,
             data = kompally_dta)
m1 %>% tidy() %>% mutate(p.value = round(p.value, 2)) 

# manually entered into regression table as lm_lin does not work with stargazer

# Table A7  --------------------------------------------------------------------
# Difference in differences models (3.6) ---

# prepare local election data
le_dta <- kompally_dta %>% 
  ungroup() %>% 
  group_by(polling_location) %>% 
  dplyr::summarise(n = sum(total_eligible),
                   n_v = sum(total_voted),
                   turnout_le = (n_v / n) * 100,
                   frt_n = sum(facial_recognition),
                   frt = ifelse(frt_n > 0, 1, 0)) 

# prepare general election data
ge_dta1 <- ge_dta %>% mutate(
  inc_vs = inc_vote_reddy / eligible_total,
  trs_vs = trs_vote_marri / eligible_total,
  bjp_vs = bjp_vote_naraparju / eligible_total,
  turnout_ge = voted_total / eligible_total ) %>% 
  select(polling_location, turnout_ge, inc_vs, trs_vs, bjp_vs)

ge_dta1 <- ge_dta1 %>% mutate(turnout_ge = turnout_ge * 100)

# create DID dataset
did_dta <- left_join(ge_dta1, le_dta, by = c("polling_location" = "polling_location")) %>% 
  dplyr::select(polling_location, frt, turnout_ge, turnout_le) %>% 
  filter(!is.na(frt)) %>% 
  pivot_longer(3:4,
               names_to = "election",
               values_to = "turnout") %>% 
  mutate(time = ifelse(election == "turnout_ge", 0, 1))

# DID model
m1 <- felm(turnout ~ frt*time | 0 | 0 | 0, data = did_dta)
m1 %>% summary(., robust = T)

stargazer(m1, 
          omit = c("Constant"),
          covariate.labels = c("Facial Recognition Technology", "Post-Period (2020 Election)", "FRT x Post"),
          se = list(m1$rse),
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          # report=('vc*s*p'),
          dep.var.labels.include = FALSE,
          keep.stat = c("n"),
          title = "Difference-in-difference effect of FRT on Turnout",
          float = F,
          star.cutoffs = c(0.1, 0.05, 0.01),
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:did",
          out = "tables/did-effect.tex") 

# Figure A7  --------------------------------------------------------------------
# Difference in differences figure (3.6) ---

did_dta <- did_dta %>% mutate(
  treatment = ifelse(frt == 1, "FRT (Treatment)", "No FRT (Control)")
)

ggplot(did_dta, aes(time, turnout, color = as.factor(treatment), linetype = as.factor(treatment))) +
  stat_summary(geom = 'line', size = 0.75) +
  scale_color_manual(values = c("black", "gray58")) +
  scale_linetype_manual(values = c("solid", "dashed")) + 
  theme_bw() +
  theme(legend.title = element_blank(), legend.position = "bottom") +
  scale_x_continuous(name = "Election",
                     breaks = c(0.00, 1),
                     labels = c("0.00" = "2019 General Election \n (Pre)", "1.00" = "2020 Local Election \n (Post)"),
                     limits = c(-0.20,1.20)) +
  scale_y_continuous(name = "Turnout")
ggsave("figures/did.png")

# Figure A8  --------------------------------------------------------------------
# Heterogenous Effects of Muslim Percent (3.7) ---

# Plot Percent Muslim
ggplot(kompally_dta, aes(x=percentmuslim_method2)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="gray58")+
  geom_density(alpha=.2, fill = "seagreen3") +
  geom_vline(aes(xintercept=median(percentmuslim_method2)),
             color="blue", linetype="dashed", size=1) +
  xlim(0,70) + 
  theme_bw() + xlab("Estimated Percent Muslim (Algorithm Approach)") + ylab("Density")
ggsave("figures/percent_muslim.png")

# Table A8  --------------------------------------------------------------------
# Interaction Effect (Continuous Effect) (3.7) ---

m1 <- felm(turnout ~ facial_recognition*percentmuslim_method2 | 0 | 0 | 0, 
           data = kompally_dta)
m1 %>% summary(., robust = T)

m2 <- felm(turnout ~ facial_recognition*percentmuslim_method2 + pfemale + total_eligible + 
             n_party + any_res | 0 | 0 | 0, 
           data = kompally_dta)
m2 %>% summary(., robust = T)

m3 <- felm(turnout ~ facial_recognition*percentmuslim_method2 | 0 | 0 | 0, 
           data = kompally_dta %>% filter(percentmuslim_method2 < 50))
m3 %>% summary(., robust = T)

m4 <- felm(turnout ~ facial_recognition*percentmuslim_method2 + pfemale + total_eligible + 
             n_party + any_res | 0 | 0 | 0, 
           data = kompally_dta %>% filter(percentmuslim_method2 < 50))
m4 %>% summary(., robust = T)

# Get exact p-values for coefficients [m1]
out <- m1 %>% summary(., robust = T)

m1_frt <- m1$coefficients[2]
m1_mus <- m1$coefficients[3]
m1_frt_mus <- m1$coefficients[4]

set.seed(1234)
mylist <- list()
for (i in 1:sims){
  mylist[[i]] <- conduct_ra(N = 36, m = 10)
}
a <- do.call("rbind", mylist) %>% as_tibble()

m1_frt_fake <- NA
m1_mus_fake <- NA
m1_frt_mus_fake <- NA

for (i in 1:sims) {
  tempdta <- a[i,] %>% t() %>% as.data.frame()
  colnames(tempdta) <- c('faketreat')
  tempdta <- cbind(kompally_dta, tempdta)
  
  mod <- lm(turnout ~ faketreat + percentmuslim_method2 + faketreat*percentmuslim_method2, data = tempdta)
  results <- coeftest(mod, vcovCL(mod, type = "HC1"))
  m1_frt_fake[i] <- results[2]
  m1_mus_fake[i] <- results[3]
  m1_frt_mus_fake[i] <- results[4]
}

# p-values for each coefficients 
m1_p1 <- mean(abs(m1_frt_fake) > abs(m1_frt)) %>% round(2)

m1_p2 <- mean(abs(m1_mus_fake) > abs(m1_mus)) %>% round(2)

m1_p3 <- mean(abs(m1_frt_mus_fake) > abs(m1_frt_mus)) %>% round(2)

# Get exact p-values for coefficients [m2]
m2_frt <- m2$coefficients[2]
m2_mus <- m2$coefficients[3]
m2_frt_mus <- m2$coefficients[8]

set.seed(1234)
mylist <- list()
for (i in 1:sims){
  mylist[[i]] <- conduct_ra(N = 36, m = 10)
}
a <- do.call("rbind", mylist) %>% as_tibble()

m2_frt_fake <- NA
m2_mus_fake <- NA
m2_frt_mus_fake <- NA

for (i in 1:sims) {
  tempdta <- a[i,] %>% t() %>% as.data.frame()
  colnames(tempdta) <- c('faketreat')
  tempdta <- cbind(kompally_dta, tempdta)
  
  mod <- lm(turnout ~ faketreat + percentmuslim_method2 + faketreat*percentmuslim_method2 +
              pfemale + total_eligible + 
              n_party + any_res, data = tempdta)
  results <- coeftest(mod, vcovCL(mod, type = "HC1"))
  m2_frt_fake[i] <- results[2]
  m2_mus_fake[i] <- results[3]
  m2_frt_mus_fake[i] <- results[8]
}

# p-values for each coefficients 
m2_p1 <- mean(abs(m2_frt_fake) > abs(m2_frt)) %>% round(2)

m2_p2 <- mean(abs(m2_mus_fake) > abs(m2_mus)) %>% round(2)

m2_p3 <- mean(abs(m2_frt_mus_fake) > abs(m2_frt_mus)) %>% round(2)

# Get exact p-values for coefficients [m3]
kd_nooutliers <- kompally_dta %>% filter(percentmuslim_method2 < 50)

m3_frt <- m3$coefficients[2]
m3_mus <- m3$coefficients[3]
m3_frt_mus <- m3$coefficients[4]

set.seed(1234)
mylist <- list()
for (i in 1:sims){
  mylist[[i]] <- conduct_ra(N = 34, m = 8)
}
a <- do.call("rbind", mylist) %>% as_tibble()

m3_frt_fake <- NA
m3_mus_fake <- NA
m3_frt_mus_fake <- NA

for (i in 1:sims) {
  tempdta <- a[i,] %>% t() %>% as.data.frame()
  colnames(tempdta) <- c('faketreat')
  tempdta <- cbind(kd_nooutliers, tempdta)
  
  mod <- lm(turnout ~ faketreat + percentmuslim_method2 + faketreat*percentmuslim_method2, data = tempdta)
  results <- coeftest(mod, vcovCL(mod, type = "HC1"))
  m3_frt_fake[i] <- results[2]
  m3_mus_fake[i] <- results[3]
  m3_frt_mus_fake[i] <- results[4]
}

# p-values for each coefficients 
m3_p1 <- mean(abs(m3_frt_fake) > abs(m3_frt)) %>% round(2)

m3_p2 <- mean(abs(m3_mus_fake) > abs(m3_mus)) %>% round(2)

m3_p3 <- mean(abs(m3_frt_mus_fake) > abs(m3_frt_mus)) %>% round(2)

# Get exact p-values for coefficients [m4]
m4_frt <- m4$coefficients[2]
m4_mus <- m4$coefficients[3]
m4_frt_mus <- m4$coefficients[8]

set.seed(1234)
mylist <- list()
for (i in 1:sims){
  mylist[[i]] <- conduct_ra(N = 34, m = 8)
}
a <- do.call("rbind", mylist) %>% as_tibble()

m4_frt_fake <- NA
m4_mus_fake <- NA
m4_frt_mus_fake <- NA

for (i in 1:sims) {
  tempdta <- a[i,] %>% t() %>% as.data.frame()
  colnames(tempdta) <- c('faketreat')
  tempdta <- cbind(kd_nooutliers, tempdta)
  
  mod <- lm(turnout ~ faketreat + percentmuslim_method2 + faketreat*percentmuslim_method2 +
              pfemale + total_eligible + 
              n_party + any_res, data = tempdta)
  results <- coeftest(mod, vcovCL(mod, type = "HC1"))
  m4_frt_fake[i] <- results[2]
  m4_mus_fake[i] <- results[3]
  m4_frt_mus_fake[i] <- results[8]
}

# p-values for each coefficients 
m4_p1 <- mean(abs(m4_frt_fake) > abs(m4_frt)) %>% round(2)

m4_p2 <- mean(abs(m4_mus_fake) > abs(m4_mus)) %>% round(2)

m4_p3 <- mean(abs(m4_frt_mus_fake) > abs(m4_frt_mus)) %>% round(2)

stargazer(m1, m2, m3, m4,
          keep = c("facial_recognition", "percentmuslim_method2", "facial_recognition:percentmuslim_method2"),
          covariate.labels = c("Facial Recognition Technology", 
                               "Percent Muslim",
                               "FRT x Percent Muslim"),
          se = list(m1$rse ,m2$rse, m3$rse, m4$rse),
          star.cutoffs = NA,
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          report=('vc*s*p'),
          dep.var.labels.include = F,
          keep.stat = c("n"),
          title = "Effect of Facial Recognition Software on Voter Turnout",
          add.lines = list(c("Mean of DV", rep(round(mean(kompally_dta$turnout), 2), 2), rep(round(mean(kd_nooutliers$turnout), 2), 4)),
                           c("Election Controls", "N", "Y", "N", "Y"), 
                           c("Polling Station Controls", "N", "Y", "N", "Y"),
                           c("Removed Outliers", "N", "N", "Y", "Y")),
          float = F,
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:het-effect-1",
          out = "tables/het-effect-1.tex") 

paste0("& [", m1_p1, "]", "$^{}$", " & [", m2_p1, "]", "$^{}$", " & [", m3_p1, "]", "$^{}$", " & [", m4_p1, "]", "$^{}$", "\\")
paste0("& [", m1_p2, "]", "$^{}$", " & [", m2_p2, "]", "$^{}$", " & [", m3_p2, "]", "$^{}$", " & [", m4_p2, "]", "$^{}$", "\\")
paste0("& [", m1_p3, "]", "$^{}$", " & [", m2_p3, "]", "$^{}$", " & [", m3_p3, "]", "$^{}$", " & [", m4_p3, "]", "$^{}$", "\\")
beepr::beep(2)

# Table A9  --------------------------------------------------------------------
# Interaction Effect (Binary Muslim Percent Variable) (3.7)) ---
kompally_dta <- kompally_dta %>% mutate(mus_med_binary = ifelse(percentmuslim_method2 > median(percentmuslim_method2), 1, 0))

m5 <- felm(turnout ~ facial_recognition*mus_med_binary | 0 | 0 | 0, 
           data = kompally_dta)
m5 %>% summary(., robust = T)

m6 <- felm(turnout ~ facial_recognition*mus_med_binary + pfemale + total_eligible + 
             n_party + any_res | 0 | 0 | 0, 
           data = kompally_dta)
m6 %>% summary(., robust = T)

m7 <- felm(turnout ~ facial_recognition*mus_med_binary | 0 | 0 | 0, 
           data = kompally_dta %>% filter(percentmuslim_method2 < 50))
m7 %>% summary(., robust = T)

m8 <- felm(turnout ~ facial_recognition*mus_med_binary + pfemale + total_eligible + 
             n_party + any_res | 0 | 0 | 0, 
           data = kompally_dta %>% filter(percentmuslim_method2 < 50))
m8 %>% summary(., robust = T)

# Get exact p-values for coefficients [m5]
m5_frt <- m5$coefficients[2]
m5_mus <- m5$coefficients[3]
m5_frt_mus <- m5$coefficients[4]

set.seed(1234)
mylist <- list()
for (i in 1:sims){
  mylist[[i]] <- conduct_ra(N = 36, m = 10)
}
a <- do.call("rbind", mylist) %>% as_tibble()

m5_frt_fake <- NA
m5_mus_fake <- NA
m5_frt_mus_fake <- NA

for (i in 1:sims) {
  tempdta <- a[i,] %>% t() %>% as.data.frame()
  colnames(tempdta) <- c('faketreat')
  tempdta <- cbind(kompally_dta, tempdta)
  
  mod <- lm(turnout ~ faketreat + mus_med_binary + faketreat*mus_med_binary, data = tempdta)
  results <- coeftest(mod, vcovCL(mod, type = "HC1"))
  m5_frt_fake[i] <- results[2]
  m5_mus_fake[i] <- results[3]
  m5_frt_mus_fake[i] <- results[4]
}

# p-values for each coefficients 
m5_p1 <- mean(abs(m5_frt_fake) > abs(m5_frt)) %>% round(2)

m5_p2 <- mean(abs(m5_mus_fake) > abs(m5_mus)) %>% round(2)

m5_p3 <- mean(abs(m5_frt_mus_fake) > abs(m5_frt_mus)) %>% round(2)

# Get exact p-values for coefficients [m6]
m6_frt <- m6$coefficients[2]
m6_mus <- m6$coefficients[3]
m6_frt_mus <- m6$coefficients[8]

set.seed(1234)
mylist <- list()
for (i in 1:sims){
  mylist[[i]] <- conduct_ra(N = 36, m = 10)
}
a <- do.call("rbind", mylist) %>% as_tibble()

m6_frt_fake <- NA
m6_mus_fake <- NA
m6_frt_mus_fake <- NA

for (i in 1:sims) {
  tempdta <- a[i,] %>% t() %>% as.data.frame()
  colnames(tempdta) <- c('faketreat')
  tempdta <- cbind(kompally_dta, tempdta)
  
  mod <- lm(turnout ~ faketreat + mus_med_binary + faketreat*mus_med_binary +
              pfemale + total_eligible + 
              n_party + any_res, data = tempdta)
  results <- coeftest(mod, vcovCL(mod, type = "HC1"))
  m6_frt_fake[i] <- results[2]
  m6_mus_fake[i] <- results[3]
  m6_frt_mus_fake[i] <- results[8]
}

# p-values for each coefficients 
m6_p1 <- mean(abs(m6_frt_fake) > abs(m6_frt)) %>% round(2)

m6_p2 <- mean(abs(m6_mus_fake) > abs(m6_mus)) %>% round(2)

m6_p3 <- mean(abs(m6_frt_mus_fake) > abs(m6_frt_mus)) %>% round(2)

# Get exact p-values for coefficients [m7]
kd_nooutliers <- kompally_dta %>% filter(percentmuslim_method2 < 50)

m7_frt <- m7$coefficients[2]
m7_mus <- m7$coefficients[3]
m7_frt_mus <- m7$coefficients[4]

set.seed(1234)
mylist <- list()
for (i in 1:sims){
  mylist[[i]] <- conduct_ra(N = 34, m = 8)
}
a <- do.call("rbind", mylist) %>% as_tibble()

m7_frt_fake <- NA
m7_mus_fake <- NA
m7_frt_mus_fake <- NA

for (i in 1:sims) {
  tempdta <- a[i,] %>% t() %>% as.data.frame()
  colnames(tempdta) <- c('faketreat')
  tempdta <- cbind(kd_nooutliers, tempdta)
  
  mod <- lm(turnout ~ faketreat + mus_med_binary + faketreat*mus_med_binary, data = tempdta)
  results <- coeftest(mod, vcovCL(mod, type = "HC1"))
  m7_frt_fake[i] <- results[2]
  m7_mus_fake[i] <- results[3]
  m7_frt_mus_fake[i] <- results[4]
}

# p-values for each coefficients 
m7_p1 <- mean(abs(m7_frt_fake) > abs(m7_frt)) %>% round(2)

m7_p2 <- mean(abs(m7_mus_fake) > abs(m7_mus)) %>% round(2)

m7_p3 <- mean(abs(m7_frt_mus_fake) > abs(m7_frt_mus)) %>% round(2)

# Get exact p-values for coefficients [m8]
m8_frt <- m8$coefficients[2]
m8_mus <- m8$coefficients[3]
m8_frt_mus <- m8$coefficients[8]

set.seed(1234)
mylist <- list()
for (i in 1:sims){
  mylist[[i]] <- conduct_ra(N = 34, m = 8)
}
a <- do.call("rbind", mylist) %>% as_tibble()

m8_frt_fake <- NA
m8_mus_fake <- NA
m8_frt_mus_fake <- NA

for (i in 1:sims) {
  tempdta <- a[i,] %>% t() %>% as.data.frame()
  colnames(tempdta) <- c('faketreat')
  tempdta <- cbind(kd_nooutliers, tempdta)
  
  mod <- lm(turnout ~ faketreat + mus_med_binary + faketreat*mus_med_binary +
              pfemale + total_eligible + 
              n_party + any_res, data = tempdta)
  results <- coeftest(mod, vcovCL(mod, type = "HC1"))
  m8_frt_fake[i] <- results[2]
  m8_mus_fake[i] <- results[3]
  m8_frt_mus_fake[i] <- results[8]
}

# p-values for each coefficients 
m8_p1 <- mean(abs(m8_frt_fake) > abs(m8_frt)) %>% round(2)

m8_p2 <- mean(abs(m8_mus_fake) > abs(m8_mus)) %>% round(2)

m8_p3 <- mean(abs(m8_frt_mus_fake) > abs(m8_frt_mus)) %>% round(2)


stargazer(m5, m6, m7, m8,
          keep = c("facial_recognition", "mus_med_binary", "facial_recognition:mus_med_binary"),
          covariate.labels = c("Facial Recognition Technology", 
                               "Above Median Percent Muslim",
                               "FRT x Above Median Percent Muslim"),
          se = list(m5$rse ,m6$rse, m7$rse, m8$rse),
          star.cutoffs = NA,
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          report=('vc*s*p'),
          dep.var.labels.include = F,
          keep.stat = c("n"),
          title = "Interaction Effect of Facial Recognition Technology and Muslim Percent on Voter Turnout",
          add.lines = list(c("Mean of DV", rep(round(mean(kompally_dta$turnout), 2), 2), rep(round(mean(kd_nooutliers$turnout), 2), 4)),
                           c("Election Controls", "N", "Y", "N", "Y"), 
                           c("Polling Station Controls", "N", "Y", "N", "Y"),
                           c("Removed Outliers", "N", "N", "Y", "Y")),
          float = F,
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:het-effect-2",
          out = "tables/het-effect-2.tex") 

paste0("& [", m5_p1, "]", "$^{*}$", " & [", m6_p1, "]", "$^{}$", " & [", m7_p1, "]", "$^{}$", " & [", m8_p1, "]", "$^{}$", "\\")
paste0("& [", m5_p2, "]", "$^{*}$", " & [", m6_p2, "]", "$^{}$", " & [", m7_p2, "]", "$^{}$", " & [", m8_p2, "]", "$^{}$", "\\")
paste0("& [", m5_p3, "]", "$^{*}$", " & [", m6_p3, "]", "$^{}$", " & [", m7_p3, "]", "$^{}$", " & [", m8_p3, "]", "$^{}$", "\\")

beepr::beep(2)

# Table A10  --------------------------------------------------------------------
# Effect of Percent Muslim on Turnout in FRT stations (3.7) ---

m1 <- felm(turnout ~ percentmuslim_method2 | 0 | 0 | 0, 
           data = kompally_dta %>% filter(facial_recognition == 1))
m1 %>% summary(., robust = T)

kompally_dta$percentmuslim_method2 %>% sd()

dv_mean <- kompally_dta %>% filter(facial_recognition == 1) %>% pull(turnout) %>% mean()

stargazer(m1,
          keep = c("percentmuslim_method2"),
          covariate.labels = c("Percent Muslim"),
          se = list(m1$rse),
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          #report=('vc*s*p'),
          dep.var.labels.include = F,
          keep.stat = c("n"),
          title = "Effect of Percent Muslim on Turnout among FRT stations",
          add.lines = list(c("Mean of DV", round(dv_mean, 2))),
          float = F,
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:turnout-percentmus",
          out = "tables/turnout-percentmus.tex") 

# Table A11  --------------------------------------------------------------------
# Effect of FRT on parties vote share (3.8) ---

m1 <- felm(BJP ~ frt_binary | 0 | 0 | 0, 
           data = ward_dta)
m1 %>% summary(., robust = T)

m2 <- felm(INC ~ frt_binary | 0 | 0 | 0, 
           data = ward_dta)
m2 %>% summary(., robust = T)

m3 <- felm(TRS ~ frt_binary | 0 | 0 | 0, 
           data = ward_dta)
m3 %>% summary(., robust = T)

# Get exact p-values for party vote share regression

set.seed(1234)
m1_p <- conduct_ri(
  formula = BJP ~ frt_binary ,
  declaration = declare_ra(N = 18, m = 6),
  assignment = "frt_binary",
  sharp_hypothesis = 0,
  data = ward_dta,
  sims = sims
)
m1_p

m2_p <- conduct_ri(
  formula = INC ~ frt_binary ,
  declaration = declare_ra(N = 18, m = 6),
  assignment = "frt_binary",
  sharp_hypothesis = 0,
  data = ward_dta,
  sims = sims
)
m2_p

m3_p <- conduct_ri(
  formula = TRS ~ frt_binary ,
  declaration = declare_ra(N = 18, m = 6),
  assignment = "frt_binary",
  sharp_hypothesis = 0,
  data = ward_dta,
  sims = sims
)
m3_p

m1_pval <- summary(m1_p)$two_tailed_p_value %>% round(2)
m2_pval <- summary(m2_p)$two_tailed_p_value %>% round(2)
m3_pval <- summary(m3_p)$two_tailed_p_value %>% round(2)

dv_bjp <- ward_dta$BJP %>% mean() %>% round(2)
dv_inc <- ward_dta$INC %>% mean() %>% round(2)
dv_trs <- ward_dta$TRS %>% mean() %>% round(2)

# Table for party effects

stargazer(m1, m2, m3, 
          keep = c("frt_binary"),
          covariate.labels = c("Facial Recognition Technology"),
          se = list(m1$rse ,m2$rse, m3$rse),
          star.cutoffs = NA,
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          report=('vc*s*p'),
          dep.var.labels.include = T,
          dep.var.labels = c("BJP Vote Share", "INC Vote Share", "TRS Vote Share"),
          keep.stat = c("n"),
          title = "Effect of Facial Recognition Software on Voter Turnout",
          add.lines = list(c("Mean of DV", dv_bjp, dv_inc, dv_trs)),
          float = F,
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:party-effect",
          out = "tables/party-effect.tex") 


paste0("& [", m1_pval, "]", "$^{}$", " & [", m2_pval, "]", "$^{}$",
       " & [", m3_pval, "]", "$^{}$", "\\")

# Table A12  --------------------------------------------------------------------
# Party difference in differences (3.8) ---

ward_ge_dta <- ge_dta %>% 
  filter(!is.na(ward_2020)) %>% 
  group_by(ward_2020) %>% 
  dplyr::summarise(inc_tot = sum(inc_vote_reddy, na.rm = T),
                   trs_tot = sum(trs_vote_marri, na.rm = T),
                   bjp_tot = sum(bjp_vote_naraparju, na.rm = T),
                   v_tot = sum(voted_total, na.rm = T)) %>% 
  mutate(INC = inc_tot / v_tot * 100,
         TRS = trs_tot / v_tot * 100 ,
         BJP = bjp_tot / v_tot * 100,
         ward = ward_2020,
         time = 0) %>% 
  select(ward, time, INC, TRS, BJP)

ge_wards <- ward_ge_dta %>% pull(ward)

ward_l_dta <- ward_dta %>% 
  filter(ward %in% ge_wards) %>% 
  mutate(time = 1) %>% 
  select(ward, time, frt = frt_binary, INC, TRS, BJP)

frt_wards <- ward_l_dta %>% filter(frt == 1) %>% pull(ward)

did_party_dta <- bind_rows(ward_ge_dta, ward_l_dta) %>% 
  mutate(frt = ifelse(ward %in% frt_wards, 1, 0))


m1 <- felm(BJP ~ frt*time | 0 | 0 | 0, data = did_party_dta)
m1 %>% summary(robust = T)

m2 <- felm(INC ~ frt*time | 0 | 0 | 0, data = did_party_dta) 
m2 %>% summary(robust = T)

m3 <- felm(TRS ~ frt*time | 0 | 0 | 0, data = did_party_dta) 
m3 %>% summary(robust = T)


dv_bjp <- did_party_dta$BJP %>% mean() %>% round(2)
dv_inc <- did_party_dta$INC %>% mean() %>% round(2)
dv_trs <- did_party_dta$TRS %>% mean() %>% round(2)

stargazer(m1, m2, m3, 
          keep = c("frt", "time", "frt:time"),
          covariate.labels = c("Facial Recognition Technology", "2020 Local Election",
                               "FRT x 2020 Local Election"),
          se = list(m1$rse ,m2$rse, m3$rse),
          #star.cutoffs = NA,
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          #report=('vc*s*p'),
          dep.var.labels.include = T,
          dep.var.labels = c("BJP Vote Share", "INC Vote Share", "TRS Vote Share"),
          keep.stat = c("n"),
          title = "Effect of Facial Recognition Software on Voter Turnout",
          add.lines = list(c("Mean of DV", dv_bjp, dv_inc, dv_trs)),
          float = F,
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:party-effect-did",
          out = "tables/party-effect-did.tex") 

# Figure A9  --------------------------------------------------------------------
# Party difference in differences plot (3.8) ---

did_party_plot <- did_party_dta %>% 
  mutate(treatment = ifelse(frt == 1, "FRT (Treatment)", "No FRT (Control)")) %>% 
  select(ward, time, treatment, INC:BJP) %>% 
  pivot_longer(4:6, names_to = "party", values_to = "vote_share")

ggplot(did_party_plot, aes(time, vote_share, 
                           color = as.factor(treatment), 
                           linetype = as.factor(treatment))) +
  stat_summary(geom = 'line', size = 0.75) +
  scale_color_manual(values = c("black", "gray58")) +
  scale_linetype_manual(values = c("solid", "dashed")) + 
  theme_bw() +
  facet_wrap(~party) + 
  theme(legend.title = element_blank(), legend.position = "bottom") +
  scale_x_continuous(name = "Election",
                     breaks = c(0.00, 1),
                     labels = c("0.00" = "2019 General Election \n (Pre)", 
                                "1.00" = "2020 Local Election \n (Post)"),
                     limits = c(-0.20,1.20)) +
  scale_y_continuous(name = "Party Vote Share",
                     limits = c(0,50))
ggsave("figures/party-vs-did.png", height = 7, width = 10)


# Table A13  --------------------------------------------------------------------
# Spillovers (3.9) ---

kompally_dta <- kompally_dta %>% 
  mutate(frt_near = case_when(polling_location == "zphs 5 kompally" ~ 0,
                              polling_location == "zphs 4 kompally" ~ 0,
                              polling_location == "zphs 3 kompally" ~ 0,
                              polling_location == "zphs 2 kompally" ~ 0,
                              polling_location == "zphs 1 kompally" ~ 0,
                              polling_location == "mpps, sakshara bharati" ~ 0,
                              TRUE ~ 1),
         frt_near_only = ifelse(frt_near == 1 & facial_recognition == 0, 1, 0)) 

m1 <- felm(turnout ~ frt_near_only | 0 | 0 | 0, 
           data = kompally_dta)
m1 %>% summary(., robust = T)

m2 <- felm(turnout ~ frt_near_only | 0 | 0 | 0, 
           data = kompally_dta %>% filter(facial_recognition == 0))
m2 %>% summary(., robust = T)

kompally_dta %>% filter(facial_recognition == 0) %>% tabyl(frt_near_only)
kompally_dta %>% tabyl(frt_near_only)

set.seed(1234)
m1_p <- conduct_ri(
  formula = turnout ~ frt_near_only,
  declaration = declare_ra(N = 36, m = 16),
  assignment = "frt_near_only",
  sharp_hypothesis = 0,
  data = kompally_dta,
  sims = sims
)
m1_p

m2_p <- conduct_ri(
  formula = turnout ~ frt_near_only,
  declaration = declare_ra(N = 26, m = 16),
  assignment = "frt_near_only",
  sharp_hypothesis = 0,
  data = kompally_dta %>% filter(facial_recognition == 0),
  sims = sims
)
m2_p

m1_pval <- summary(m1_p)$two_tailed_p_value %>% round(2)
m2_pval <- summary(m2_p)$two_tailed_p_value %>% round(2)

dvmean1 <- kompally_dta %>% pull(turnout) %>% mean() %>% round(2)
dvmean2 <- kompally_dta %>% filter(facial_recognition == 0) %>% pull(turnout) %>% mean() %>% round(2)

stargazer(m1, m2, 
          keep = c("frt_near_only"),
          covariate.labels = c("Near FRT Station"),
          se = list(m1$rse, m2$rse),
          star.cutoffs = NA,
          digits = 2,
          # ci = T,
          # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
          report=('vc*s*p'),
          dep.var.labels.include = F,
          dep.var.caption = "DV: Turnout",
          keep.stat = c("n"),
          title = "FRT Spillovers",
          add.lines = list(c("Mean of DV", dvmean1, dvmean2)),
          float = F,
          table.placement = "H",
          omit.table.layout = "n",
          notes.append = F,
          no.space = T,
          label = "table:spillovers",
          out = "tables/spillover-effect.tex") 

# Replace regular p-values with exact p-values (Fisher)
paste0("& [", m1_pval, "]", "$^{}$", " & [", m2_pval, "]", "$^{}$", "\\")

